----------general least squares fit------
-- finds coefficients a[i] for best fit of weighted sum of functions of x
include "::incl:matrix ops"

fit(x) = sum(a[i,1]*f(x)[i],i,1,nparms)

f(x)[j]=x^(j-1)  -- polynomial a1+a2*x+a3*x^2...

P[i,j]=f(x[i])[j] dim[ndata,nparms]
PtP:=multiply(transpose(P),P):
Pty:=multiply(transpose(P),y):
inv:=invert(PtP):
a:=multiply(inv,Pty):

data=read(xydata)
x[i]=data[i,1]; y[i]=data[i,2] dim[ndata]; 

ndata=count(data)
nparms=3

a:{{5.582},{0.159},{0.006}}

plot data
plot fit(X)
